The goal of this notebook is to phase more clusters of variants using
the same approach we used to define the ‘genotypes’
Genome-1 and Genome-2 in the previous
notebook. The output of this notebook is the list of variants
annotated by whether they are part of any
haplotypes/clusters.
The inputs for this analysis are (1) the filtered and combined
variants annotated by major genotypes from
02-determine-main-genotypes.Rmd and (2) annotations of the
positions of the main MeV open reading frames.
## ==== File paths input data ==== ##
if (exists("snakemake")) {
# Variant data annotated with major genotypes
labeled.data = snakemake@params[['incsv']]
# Annotations for the standard Measles ORFs
annotations.filepath = snakemake@params[['annotations']]
# Path to the ancestral mutations
ancestral.data = snakemake@params[['ancestral']]
} else {
# Variant data annotated with major genotypes
labeled.data = "../../results/variants/genotyped_variants.csv"
# Annotations for the standard Measles ORFs
annotations.filepath = "../../config/annotations.csv"
# Path to the ancestral mutations
ancestral.data = "../../config/ref/annotated_SSPE_consensus_snps.csv"
}
## ==== File paths output data ==== ##
# Path to save the results
if (exists("snakemake")) {
# Variants annotated by new haplotype clusters
output.path = snakemake@params[['outcsv']]
# Path to save figures
figure.dir = paste0(snakemake@params[['figures']], "/")
} else {
# Variants annotated by new haplotype clusters
output.path = "../../results/variants/clustered_variants.csv"
# Path to save figures
figure.dir = "../../results/figures/"
}
# Make the figure directory if it doesn't exist
if (!file.exists(figure.dir)) {
dir.create(figure.dir, recursive = TRUE)
print("Directory created.")
} else {
print("Directory already exists.")
}
## [1] "Directory already exists."
Since the goal of this notebook is use the same technique that I used
to identify Genome-1, Genome-2, and
Genome-1-1, but for mutations that aren’t necessarily
found in every single tissue, I need to account for missing
mutations. To do this, I’ll spread the data frame of frequencies and
assign an allele frequency of 0 if a variant is missing from a tissue.
This will allow us to cluster all of the mutations, even if they weren’t
observed in a tissue sample.
I’ll write the approach I used in
02-determine-main-genotypes.Rmd into a reusable function. I
can provide a list of SNPs, the full data frame, and the number of
clusters I’m targeting.
cluster.snps <- function(list.of.snps, snp.df, n.clusters) {
# SNP as column and Allele frequency as row
frequency.by.snp = snp.df %>%
filter(SNP %in% list.of.snps) %>%
select(AF, SNP, Tissue) %>%
pivot_wider(names_from = SNP, values_from = AF, values_fill = 0) %>%
select(!Tissue)
# Calculate R between every pair of columns while handling NAs
snp.correlation = cor(frequency.by.snp, use="pairwise.complete.obs")
# Convert to distance (positive corr is close to 0)
snp.dist = as.dist(1 - snp.correlation)
# Create k-medoids clustering with n clusters
snp.kmedoids = pam(snp.dist, n.clusters)
kclusters = snp.kmedoids$cluster
# Convert to a data frame
kclusters.df = data.frame(SNP = names(kclusters), cluster = kclusters)
# Assign the clusters to the original data frame
kmedoids.SNPs = snp.df %>%
filter(SNP %in% list.of.snps) %>%
left_join(., kclusters.df, by = "SNP") %>%
mutate(cluster = if_else(is.na(cluster), "no cluster", as.character(cluster)))
# Sort the cluster names by mean frequency
cluster.order = kmedoids.SNPs %>%
group_by(cluster) %>%
summarize(AF = mean(AF)) %>%
arrange(-AF) %>%
mutate(order = 1:n.clusters) %>%
select(cluster, order)
kmedoids.SNPs = kmedoids.SNPs %>%
left_join(., cluster.order, by = "cluster") %>%
select(!cluster) %>%
rename(cluster = order) %>%
mutate(cluster = as.character(cluster))
# Add the number of SNPs per cluster
n.clusters.per.snp = kmedoids.SNPs %>%
select(SNP, cluster) %>%
distinct() %>%
group_by(cluster) %>%
count() %>%
mutate(cluster_size = paste0(cluster, " (", n, " snps in cluster)"))
return(left_join(kmedoids.SNPs, n.clusters.per.snp, by = "cluster"))
}
Although this approach works reasonably well, it’s very susceptible to noise. This means that you can’t just provide all the variants at once and expect the algorithm to partition them into correct clusters. Instead, the best approach I found was to break variants down into bins based on their frequency and cluster these separately. The higher average frequency, the less noise, and the better the algorithm is at partitioning the variants into meaningful clusters.
First, I’ll start with a ‘high-frequency’ bin. This includes all variants that reach greater than or equal to 25% in at least one tissue.
Some clusters have only a single mutation in them. Clearly, some of
these belong to Genome-1, Genome-2, and
Genome-1-1, but are missing from a single tissue, perhaps
due to low coverage in that particular tissue.
Cluster 5 (L-A15084G) is a mutation in Genome-1-1 that
was too low frequency to be called in the Frontal Cortex 2
sample. Cluster 6 (M-C4502T) is mutation that’s part of
Genome-2 but was too low coverage in the
Cerebellum. Finally, Cluster 9 (M-C4573T) is a mutation in
Genome-1 that was also too low coverage in the
Cerebellum. Generally, the Cerebellum has the
lowest coverage of any tissue, so this isn’t surprising. All-in-all,
this adds three mutations to already existing haplotypes. I already
addressed F-C7036T, N-G810A, and H-T7293C in the
identify-backgrounds.Rmd notebook. Finally, M-C4532T might
also belong as part of the reference?
I’ll annotate these variants in the main data frame.
# Clusters that are clearly genome-1
genome.1.missing.snps = twentyfive.percent.or.more.clusters.df %>%
filter(cluster %in% c(4)) %>%
pull(SNP) %>%
unique(.)
# Clusters that are clearly genome-1-1
genome.1.1.missing.snps = twentyfive.percent.or.more.clusters.df %>%
filter(cluster %in% c(5)) %>%
pull(SNP) %>%
unique(.)
# Clusters that are clearly genome-2
genome.2.missing.snps = twentyfive.percent.or.more.clusters.df %>%
filter(cluster %in% c(6)) %>%
pull(SNP) %>%
unique(.)
# Clusters that break the rules
maybe.in.both.snps = twentyfive.percent.or.more.clusters.df %>%
filter(cluster %in% c(3, 7, 2, 1)) %>%
pull(SNP) %>%
unique(.)
# Annotate the labeled data frame with what's possible here
updated.haplotype.df = labeled.df %>%
mutate(Background = case_when(SNP %in% genome.1.missing.snps ~ "genome-1",
SNP %in% genome.1.1.missing.snps ~ "genome-1",
SNP %in% genome.2.missing.snps ~ "genome-2",
SNP %in% maybe.in.both.snps ~ "both",
TRUE ~ Background)) %>%
mutate(Haplotype = case_when(SNP %in% genome.1.missing.snps ~ "genome-1",
SNP %in% genome.1.1.missing.snps ~ "genome-1-1",
SNP %in% genome.2.missing.snps ~ "genome-2",
SNP %in% maybe.in.both.snps ~ "both",
TRUE ~ Haplotype))
Now, lets look at the remainder of the clusters of SNPs. These are SNPs that actually clustered with other variants.
All of these look like reasonable clusters, so I’ll annotate the
variant data frame with these new clusters (subclonal haplotypes). I’ll
also try to guess if these are on the background of
Genome-1 or Genome-2 if it’s possible.
cluster.1.snps = twentyfive.percent.or.more.clusters.df %>%
filter(cluster %in% c(8)) %>%
pull(SNP) %>%
unique(.)
# The missing SNPs in Cerebellum are due to coverage
cluster.2.snps = twentyfive.percent.or.more.clusters.df %>%
filter(cluster %in% c(11, 14)) %>%
pull(SNP) %>%
unique(.)
cluster.3.snps = twentyfive.percent.or.more.clusters.df %>%
filter(cluster %in% c(10)) %>%
pull(SNP) %>%
unique(.)
cluster.4.snps = twentyfive.percent.or.more.clusters.df %>%
filter(cluster %in% c(13)) %>%
pull(SNP) %>%
unique(.)
cluster.5.snps = twentyfive.percent.or.more.clusters.df %>%
filter(cluster %in% c(15)) %>%
pull(SNP) %>%
unique(.)
cluster.6.snps = twentyfive.percent.or.more.clusters.df %>%
filter(cluster %in% c(16)) %>%
pull(SNP) %>%
unique(.)
cluster.7.snps = twentyfive.percent.or.more.clusters.df %>%
filter(cluster %in% c(17)) %>%
pull(SNP) %>%
unique(.)
cluster.8.snps = twentyfive.percent.or.more.clusters.df %>%
filter(cluster %in% c(18)) %>%
pull(SNP) %>%
unique(.)
cluster.9.snps = twentyfive.percent.or.more.clusters.df %>%
filter(cluster %in% c(19)) %>%
pull(SNP) %>%
unique(.)
cluster.10.snps = twentyfive.percent.or.more.clusters.df %>%
filter(cluster %in% c(20)) %>%
pull(SNP) %>%
unique(.)
# Annotate the labeled data frame with what's possible here
updated.haplotype.df = updated.haplotype.df %>%
mutate(Haplotype = case_when(SNP %in% cluster.1.snps ~ "cluster 1",
SNP %in% cluster.2.snps ~"cluster 2",
SNP %in% cluster.3.snps ~"cluster 3",
SNP %in% cluster.4.snps ~ "cluster 4",
SNP %in% cluster.5.snps ~ "cluster 5",
SNP %in% cluster.6.snps ~ "cluster 6",
SNP %in% cluster.7.snps ~ "cluster 7",
SNP %in% cluster.8.snps ~ "cluster 8",
SNP %in% cluster.9.snps ~ "cluster 9",
SNP %in% cluster.10.snps ~ "cluster 10",
TRUE ~ Haplotype))
Now, I’ll move onto the next bin. These are mutations that reach at least 5% in one tissue but weren’t part of the previous analysis. Most of these will be challenging to haplotype due to noise. I’ll only try to specify the clusters that seem the most clear cut.
Most of these aren’t very obviously clusters that represent true haplotypes. They’re more likely just an artifact of noise. I took the 8 that looked the most promising. It’s probably possible to do this more systematically using the internal variation in a cluster, but the results will be the same and there aren’t enough clusters to justify the need for this approach.
cluster.11.snps = twentyfive.percent.or.fewer.clusters.df %>%
filter(cluster %in% c(15)) %>%
pull(SNP) %>%
unique(.)
cluster.12.snps = twentyfive.percent.or.fewer.clusters.df %>%
filter(cluster %in% c(4)) %>%
pull(SNP) %>%
unique(.)
cluster.13.snps = twentyfive.percent.or.fewer.clusters.df %>%
filter(cluster %in% c(14)) %>%
pull(SNP) %>%
unique(.)
cluster.14.snps = twentyfive.percent.or.fewer.clusters.df %>%
filter(cluster %in% c(12)) %>%
pull(SNP) %>%
unique(.)
# Annotate the labeled data frame with what's possible here
updated.haplotype.df = updated.haplotype.df %>%
mutate(Haplotype = case_when(SNP %in% cluster.11.snps ~ "cluster 11",
SNP %in% cluster.12.snps ~ "cluster 12",
SNP %in% cluster.13.snps ~ "cluster 13",
SNP %in% cluster.14.snps ~ "cluster 14",
TRUE ~ Haplotype))
We were able to cluster a reasonable number of the high-frequency mutations (>= 5% frequency). How many mutations haven’t been haplotyped yet?
## [1] "There are about 471 that are still totally un-haplotyped."
## [1] "There are 60 that have been clustered."
## [1] "There are still about 262 that can reasonably be haplotyped (Above 5% in more at least one sample)."
Here are the frequency of all current haplotypes in the 15 tissue samples.
Looking closely at the haplotypes, there is one amendment to be made. Cluster 10 and 14 are likely the same cluster. I think these weren’t clustered together because the mean frequency of these mutations is split between the frequency bins that I used to partition mutations.
For now, I’ll combine these into a single cluster. I’ll explore this further with bridging reads later in the analysis.
updated.haplotype.df = updated.haplotype.df %>%
mutate(Haplotype = case_when(Haplotype %in% c("cluster 14") ~ "cluster 10",
TRUE ~ Haplotype))
I’ll write out these newly phased clusters of mutations for further analysis.
## [1] "Writing the results to results/variants/clustered_variants.csv"
At this stage in the analysis, we’ve phased mutations into rough
haplotypes. We’ve also assigned some mutations to the major ‘genotypes’
(Genome 1 and Genome 2) that were missed in
the previous notebook due to low coverage.
Let’s visualize all of these mutations in each tissue specimen.
Mutations are annotated by whether they’re Genome 1,
Genome 2, or a ‘Driver’ mutations. These ‘Drivers’ are
mutations with impacts similar to mutations that have been observed in
previous cases of SSPE.
What are the SNPs we’re calling ‘Genome 1 (G1)’. These will be different than the ‘Candidate Genome 1’ SNPs in Figure 2, because these were identified by their correlation in frequency between all tissue specimens.
In the previous notebook, we noticed a set of mutations that made up
a haplotype that arises on the background of Genome 1, but
it mostly missing from the Frontal Cortex 2 sample. We called this
haplotype Genome-1-1, and it’s very significant because
it’s indicative some early divergence in the frontal cortex. Let’s take
a look at this haplotype in every tissue.